%%  mf_build
% All the cells of the finite difference grid are given properties.

% Todo:
% - !injected water in input doesnt seem to correspond with output
% - remove all rech-lines
% - convert to real size
% - re-add rech package
% - apply trans-evap package

clear variables; close all; clc;

%%  Set variables
% load CurrentSeriesNumber n
load curNum;    % current setup number: n
load setup setup; % load setup struct

freshwater  = 0; % density difference compared to water; [kg/m3]
saltwater   = 1; % density difference compared to water; [kg/m3]

ixWall  = setup(n).grExtra.ixWall;  % Convert to IBOUND, so nodegrid is used
ixRech  = setup(n).grExtra.ixRech;  % Ceil/floor consturtion is used, so RECH is always on the leftside of the wall

%% Model arrays located on the picture, cell properties are being added through IBOUND matrix to the cellnodes
IBOUND  = setup(n).gr.const(2);  % Define active cells
K       = setup(n).setupVars.CONDUC/(24*60);  % CONDUC[cm/d]->[cm/min]
HK      = setup(n).gr.const(K);
VK      = HK;
VK(:,setup(n).grExtra.ixhInjB,2)=0;
PEFF    = setup(n).gr.const(setup(n).setupVars.POREFF);
PEFF(:,1:2,1:2)=1;

MH      = max(setup(n).gr.zGr);
STRTHD  = setup(n).gr.const(0);

%%  Set initial concentration
ICBUND = setup(n).gr.const(saltwater);
ICBUND(IBOUND==-1)=-1;

if setup(n).setupVars.REUSECONC == 1            % if true the final concentration of the previous cycle is beeing used
    STCONC = setup(n-1).results.CsavedTresh;
else
    STCONC = setup(n).gr.const(saltwater);      % start concentration is all saltwater
end

%%  Constand Head Package
CHDzone             = 3;
IBOUND(:,setup(n).gr.Nx,end) = CHDzone;  % CHD package
CHDDENSOPT          = 2;        % env. hd at ocean boundary, Langevin et al 2008, p22
[CHD, PNTSRC]       = setup(n).gr.bcnZone(setup(n).basename,'CHD',IBOUND,[CHDzone  0 0 CHDDENSOPT],saltwater);

%%  RECH [L/T]
% The injected water is named milk
if 0
    IBOUND(:,1:ixRech-1,1) = 6;                   % RECH package
    RECH    = zeros(setup(n).gr.Ny,setup(n).gr.Nx,length(setup(n).perdata.perlen));
    [NrMilk_y NrMilk_x NrMilk_z] = find(IBOUND==6);
    RECHarea= sum(setup(n).gr.AREA(NrMilk_x));
    QMilk   = setup(n).perdata.milkSign*setup(n).setupVars.MULTIMILK;      % Total amount to be injected [cm3/min]
    
    for a = 1:length(setup(n).perdata.perlen);
        for b = 1:length(NrMilk_x)              % For all rech coordinates
            RECH(NrMilk_y(b),NrMilk_x(b),a) = QMilk(a).*(1/RECHarea); % since just a part is modeled
        end
    end
    setup(n).setupVars.RECH   = RECH;
end

%%  Wells
[well,WEL,PNTSRC(end+1:end+length(setup(n).perdata.perlen)),NPER]=mf_setwells(setup(n).basename,setup(n).gr,HK,'wells');
setup(n).well   = well;

%%  Horizontal Flow Barrier
% ILay    = 1:abs(setup(n).gr.zm);  % Layer index
ILay    = 1:setup(n).grExtra.izWall;
HFBc    = 10000;            % the hydraulic resistance per layer, only one layer is used in zoetwatervat
HFB1    = mf_setHFB([setup(n).setupVars.XWALL,setup(n).gr.ym(1);setup(n).setupVars.XWALL,setup(n).gr.ym(end)],setup(n).gr.xGr,setup(n).gr.yGr,ILay,HFBc);
HFB2    = mf_setHFB([setup(n).setupVars.XWALL,setup(n).gr.ym(1);setup(n).setupVars.XWALL,setup(n).gr.ym(end)],setup(n).gr.xGr,setup(n).gr.ym,ILay,HFBc);  % Some inconsistency within interpolating the y-values

ILay2   = setup(n).grExtra.izvInjB;
% injectionbound
HFB3    = mf_setHFB([setup(n).gr.xm(setup(n).grExtra.ixvInjB),min(setup(n).gr.ym(setup(n).grExtra.iyvInjB));setup(n).gr.xm(setup(n).grExtra.ixvInjB),max(setup(n).gr.ym(setup(n).grExtra.iyvInjB))],setup(n).gr.xGr,setup(n).gr.yGr,ILay2,HFBc);     % Vertical injection boundary
HFB4    = mf_setHFB([setup(n).gr.xm(setup(n).grExtra.ixvInjB),min(setup(n).gr.yGr(setup(n).grExtra.iyvInjB));setup(n).gr.xm(setup(n).grExtra.ixvInjB),max(setup(n).gr.yGr(setup(n).grExtra.iyvInjB))],setup(n).gr.xGr,setup(n).gr.ym,ILay2,HFBc);     % Vertical injection boundary

HFB     = [HFB1;HFB3;HFB2;HFB4];      % Needed for mf_setup
setup(n).HFB = HFB;         % Needed for mf_analyze and archiving

%%  Save variables
save setup setup
if isfield(setup,'outputPath'),  save(strcat(setup(n).outputPath,'setup'),'setup'); end
%%  Save variables to keep compatibility with mf_setup
basename = setup(n).basename;       % Needed for mf_setup
gr       = setup(n).gr;                % Needed for mf_setup